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Global dynamical behaviors of the competitive Lotka-Volterra system even in 3-dimension are not 
fully understood. The Lyapunov function can provide us such knowledge once it is constructed. In 
this paper, we construct explicitly the Lyapunov function in three examples of the competitive Lotka- 
Volterra system for the whole state space: (1) the general 2-dimensional case; (2) a 3-dimensional 
model; (3) the model of May-Leonard. The dynamics of these examples include bistable case 
and cyclical behavior. The first two examples are the generalized gradient system defined in the 
Appendixes, while the model of May-Leonard is not. Our method is helpful to understand the limit 
cycle problems in general 3-dimcnsiorial case. 
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I. INTRODUCTION 

Lotka-Voltcrra system is one of the most fundamental 
models describing the interaction of n species in mathe- 
matical ecology [1] , physics and economics [2] . It is given 
by the following ordinary differential equations: 

Definition 1 (Lotka-Volterra system). 

Xi = Xi - ^ aijXj^ ,i = l,...,n, (1) 

where each Xi{i = 1, • • • ,n) represents the population of 
one species and bi,aij{i = 1, • • • ,n;j = I,-- - ,n) are 
constants depending on the environment. The state space 
of the system (1) is represented by the non-negative vec- 
tors R\ = {{xi,--- ,Xn) e > 0, i = 1, • • • , n}. 
When hi > 0, aij > 0(i = 1, • • • , n), it is the competitive 
Lotka- Volterra system. 

Due to the nonlinear attributes of the competitive 
Lotka-Volterra system, its dynamics can be complex 
when n > 3, such as cyclical behavior [3] and chaotic 
behavior [4]. M. Hirsch has proved any trajectory of a n- 
dimensional competitive Lotka-Volterra system will con- 
verge to an invariant surface E, homeomorphic to (n — 1)- 
dimcnsional unit simplex A = Xj : Xj > 0, 5^"^/ x.^ = 1. 

In 3-dimensional case, following M. Hirsch's general 
result, M. L. Zeeman identified 33 stable equivalence 
classes, of which only classes 26 — 31 can have limit cy- 
cles. Then in [5], D. Xiao and W. Li proved the number 
of limit cycles is finite without a heteroclinic polycycle. 
In [6], J. Hofbaucr and J. So conjectured the number of 
limit cycles is at most two. Nevertheless, three limit cy- 
cles were constructed numerically in [7, 8] and four in 
[9, 10]. M. L. Zeeman also tried to deduce global dy- 
namics by the edges of the carrying simplex [11, 12]. Till 



now, however, the question of how many limit cycles can 
appears in M. L. Zeeman's six classes 26 — 31 remains 
open. 

To analyze global dynamics, M. Planck studied the 
Lotka-Volterra system by hamiltonian theory, however, 
in limited parameter region [13]. The split Lyapunov 
function has been used [14] , but it is not monotone along 
all the trajectories in the state space, hence constructed 
locally. Besides, the classical Lyapunov function has also 
been constructed for n-dimensional case in the parameter 
region with one global stable equilibrium [15, 16]. As 
there is no general way of constructing the Lyapunov 
function [17], this theory has not been explored more to 
study the competitive Lotka-Volterra system. 

In this paper, based on the framework of general dy- 
namics recently proposed in [18, 19], we construct the 
Lyapunov function in three examples of the compet- 
itive Lotka-Volterra system. This Lyapunov hmction 
is monotone along trajectories in the state space, thus 
can demonstrate global dynamics. The first example is 
the general 2-dimensional case [15]. The second is a 3- 
dimensional system given by [15] and the construction 
method is the same with the first's as both are the gen- 
eralized gradient system defined in the Appendixes. The 
third example is the classical May-Leonard 3-dimensional 
system [3]. As it is not the generalized gradient system, 
we provide there a different construction method. 

This paper is organized as follows. In Sec. II, we 
uniformly construct the Lyapunov function for the 2- 
dimensional model, and analyze its dynamics in the state 
space. In Sec. Ill, we study the 3-dimensional competi- 
tive Lotka-Volterra system given by [15]. In Sec. IV, we 
study the model of May-Leonard. In Sec. V, we summa- 
rize our work. In the Appendixes, we introduce briefly 
our construction framework, discuss the generalized gra- 
dient system, and then give detailed calculation on other 
dynamical parts in our framework of the three examples. 
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II. THE GENERAL 2-DIMENSIONAL 
COMPETITIVE LOTKA-VOLTERRA SYSTEM 

Example 1. The general 2 -dimensional competitive 
Lotka-Volterra system is given by: 

Xi = Xi{bi — Xi — ax2) 
X2 = X2{b2 - f3xi - X2) ' 



(2) 



where bi, 62, a, j3 are non-negative constants [15]. 
By setting Xi = X2 = 0, four non-negative equi- 
libriums are derived: (1) a positive one = 
{b\ — ab2,b2 — I3bi) /(I — a/3) existing when a < 61/62, 
13 < 62/61 ora> 61/62, /3 > 62/61; (2) E+o = (61, 0); 
(3) Eo+ = (0,62); (4) Eoo = (0,0). Here the subscript + 
denotes the population of the species is positive and the 
subscript means the species dies out. 



A. Construction of the Lyapunov function 

Now we introduce our method to construct the Lya- 
punov function of the system. This method is based on 
the framework in [18, 19]. The idea is as follows. Assume 
there is a Lyapunov function (f) and its partial derivative 
is given by 



Aii{xi,X2) Ai2{xi,X2) \ f Xi 
A2l{xi,X2) A22{X1,X2) ) \ ±2 



where Aii{xi,X2), ^12(2:1, X2), ^2i(-i;i, -^2), ^22(2:1,. T2) 
are undetermined coefficients. Our aim is to choose 
proper coefficients so that: (1)V x V(p = 0; (2)^ < 0, 
i.e., Lie derivative of (j) decreasing along trajectories. 

We discover that Aii(a;i,a;2) = (3/xi, Ai2{xi,X2) = 0, 
^21(2:1, 0:2) = 0, A22{xi,X2) = a/x2 is a proper setting. 
Thus we get 



(3) 



1 ^ = -«(^2-M-X2) ■ 

With direct calculation, V x = and 

; d(j) . dcj) . 
9 = ^—xi + ^—X2 
0x1 0X2 

= -Pxi{bi - Xi - 0X2)^ - (XX2{f)2 — Pxi - X2)^ < , 

as x\ and X2 are all non-negative population species and 

f3 and a. are all non-negative constants, (/'(x) = hap- 
pens only at X e UsgR2^w(s), where ix)(s) denotes the 
w-limit set [20]. Thus, we can get a Lyapunov function 
by integrating the Eq. (3): 

Q Q, 

(\) = '-t\ + -X2 - (3biXi - 062X2 -I- a^a;iX2 . (4) 

Here we mention that the choice on the coefficients 
Aix{,xi,X2), ^i2(a;i,a;2), A2\{xi,X2), ^22(2:1, 0:2) is not 
unique. Our choice is straightforward and meets the re- 
quirements. 



a^(l - ap) . (5) 



B. Analysis on dynamics in the state space 

In this subsection, we give a classified discussion on 

dynamics in each parameter region by the Hessian matrix 
of the Lyapunov function at As = /3, = 

^' dxidx2 ~ dx2dxi ~ ^® ^^"^ determinant of 
Hessian matrix: 

dx\ 8x2 \dx1dx2. 

Thus, the type of dynamics can be classified into four 
cases: 

(1) Stable coexistence case: a < 61/62, /3 < 62/61. 

A > and > indicate is a globally stable 
equilibrium with the minimum energy value. 

(2) Bistable case: a > 61/62, 13 > 62/61. 

A < indicates £■++ is a saddle point. As the system 
(2) is bounded in the first quadrant, it has two stable 
equilibriums £^+0 and £'0+ on the boundary. 

(3) One survival case: a < 61/62, /3 > 62/61 or a > 61/62, 

/?< 62/61. 

It has one globally stable equilibrium on an axis of 
coordinate, E^q appears when a < 61/62, /3 > 62/61 
or i?o+ appears when a > 61/62, (3 < 62/61. We just 
show the case where the species xi survives in Fig. 1, 
i.e., when a < 61/62, /3 > 62/61. The case where the 
species X2 survives can be shown similarly. 

(4) Degenerate case: a = 61/62, /3 = 62/61. 

The Lyapunov function has the minimum value along 
the line: 

\/6i62 — 1/62/61X1 — 1/61/62X2 = as in this case 



4'=^ (^Vhh - Vh/bixi - ^61/62X2^ 



- 6162 . (6) 



Each trajectory will converge to one of the points on 
the line, depending on the initial value. 

Four remarks arc made here: 

• Our result on dynamics of the system is consistent 
with the stability analysis near equilibriums in [15]. 

Additionally, our Lyapunov function is constructed 
uniformly for the whole parameter space and thus 
can provide dynamics for any perturbation on the 
parameters. Therefore, the criterion on the classi- 
fication on the type of dynamics due to parameters 
changing can be based on the Lyapunov function: 
(1) when exists, the Hessian matrix of the 

Lyapunov function at £++ can show its type of 
stability and we have the previous two cases; (2) 
when does not exist, we have the case (3); (3) 
the remaining one is the case (4). 
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Case (3) Case (4) 



FIG. 1. The energy landscape of example (1) for the various 
cases: (1) a = /3 = ^, bi — b2 — 1: stable coexistence case; 
(2) Q = ^ = 2, &i = 62 = 1: bistable case; (3) q = /3 = 2, 
bi — h2 — 1: one survival case; (4) a = /3 = fei = 62 = 1: 
degenerate case. 

• Visualization with the energy landscape of the Lya- 
punov function (Fig. 1) provides a clear observation 
on dynamics in each case above: Fig. I's case (1) 
to case (4) relate to these cases (1) to (4) discussed 
above respectively. Saddle-node bifurcation can be 
observed from this energy landscape. The bifurca- 
tion happens from the case (1) to the case (2) and 
the degenerate case (4) has the minimum energy 
value on a line. 

• We show here the dynamics of the system can di- 
rectly be determined by the Lyapunov function 
alone. In [21], M. L. Zeeman used the Lyapunov 
function to prove that the stable nullcline classes 
coincide with the stable topological classes in this 
system. 

• In the Appendixes, we will define the general- 
ized gradient system, as a natural generalization 
to typical gradient system. We thus find this 2- 
dimensional system meet the definition. Further- 
more, the construction method can be applied to 
the generalized gradient system in n-dimension, 
and we will show a 3-dimensional example in next 
section. 



III. A 3-DIMENSIONAL MODEL 

Following the definition of the generalized gradient sys- 
tem in the Appendixes, we find that a 3-dimensional com- 



petitive Lotka-Volterra system given by [15] is another 
one and thus we construct the Lyapunov function by the 
same method as the last section's. 

Example 2. 

Xi = Xi (1 — xi — ax2) 
±2 = 2:2 (1 - (ixi - X2 - I3xz) , (7) 
±3 = 2^3 (1 - ctX2 - 2:3) 

where a,/3 are non-negative coefficients. By setting Xi — 
X2 — X3 = 0, non-negative equilibriums are derived: (1) 
a positive one = (1 — q, 1 — 2/?, 1 — a) /(I — 2a/3) 

existing when (1 — a) /{I — 2a/3) > and (1 — 2/3)/(l — 
2a(i) > 0; (2) E+o+ = (1,0, 1); (3) Eo+o = (0, 1,0); (4) 

i;ooo = (0,0,0). 

A. Construction of the Lyapunov function 

With the similar method used in the general 2- 
dimensional competitive Lotka-Volterra system, we 
choose the corresponding undetermined matrix to be 

l3/xi \ 
a/x2 
f3/x3 J 

then 

' ||^_^(l_^^_aX2) 

< ^ = -a{l-l3xi-X2-f3x3) . (8) 
, si = -/3(1 -axa -X3) 

Since V x Wcj) ~ and the Lie derivative of (j> is 

4> — — Pxi (1 — xi — 0x2)^ — ax2 (1 — jixi — X2 — Px^y^ 
- /3x3 (1 - ax2 - xsf < , 

as a;i,a;2 and x^ are all non- negative population species 
and (3 and a are all non-negative constants. 0(x) = 
happens only at x G UggRS w(s). We can construct a 
Lyapunov function by integrating the Eq. (8). 

^ix'i+xl)-\-^xl-\-a(3{xiX2+X2X3)~(3{xi+X3)-ax2 ■ 

(9) 

B. Analysis on dynamics in the state space 

With the Lyapunov function constructed globally on 
R^, first the classified stability analysis near equilibriums 
by [15] can be unified now. Second, when a = 1 and 
/3 = l/2, 

= J [{xi +X2- if + (X2 -f - 1)' - 2] (10) 
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indicates the degenerate case of the system has the mini- 
mum energy on the intersection of the surfaces X1+X2 — 
1 = and X2+XS — I = 0. Third, as it is a generaUzed gra- 
dient system without a trajectory contouring along the 
energy landscape of the Lyapunov function, the system 
(7) does not have a limit cycle [19] . 



IV. THE MODEL OF MAY-LEONARD 

In [3], R. May and W. Leonard studied a 3-dimensional 
competitive Lotka-Volterra system. 



Example 3. 



ii = xi {1 — xi — ax2 — 13x3) 
X2 = 2:2 (1 - I3xi — X2~ ax^) 
X3 = (1 - axi - I3x2 - X3) 



(11) 



where a, (3 are non-negative coefficients. The possi- 
ble equilibriums contain: (1)(0,0,0); (2) three single- 
population survive (1,0, 0), (0, 1,0), (0, 0, 1); (3) three 
two-population solutions of the form ^^""^"'^^'"^ ; (4) and 
three-population survive 



(1.1,1) 

l+a+/3' 



A. Construction of the Lyapunov function 

We find that this model is not the generalized gradient 

system, and thus the construction method in Sec. II can 
not be applied here. So we will give another method to 
construct the Lyapunov function in this section. 

For convenience, let us introduce some new variables: 
7 = a + /3 — 2, P — X1X2X3 and O = Xi -\- X2 + x^. Then 

P = ±1X2X3 + XlX2X:i + X1X2X3 

= P [3 - (1 + a + /3)0] 
= P[3-(3 + 7)0] 
= P[3(l-0)-70], 

and 

O = Xx-\- ±2+ Xz 

= - [xl -\- xl -\- xl -\- {a + j3){xiX2 + X2X3 + X3X1)] 

= 0(1 - O) - j{xiX2 + X2X3 + X3X1), 

where the P and O denotes the Lie derivative of P and O 
respectively. Next, we construct the Lyapunov function 
in two different parameter regions: (1)7 = 0; (2)7 ^ 0. 

1. When 7 = 0: 

Noting that P = 3P(1 - O) and O = 0(1 - O), 
thus 



--+30 = -3(1 -0)2 <0. 



(12) 



So if we can a function whose Lie derivative is 
—P/P + 30, then it is a Lyapunov function. This 
can be done by simply integrate — P/P + 30 and 
we get a Lyapunov function 

= 3O-lnP (13) 
= 3(a;i + a;2 + 0:3) - ln(a;ia;2a;3). (14) 

(p = only when 0—1 = 0, i.e., all the trajectories 
converge to the plane = 1. 

2. When 7^0: 

Noting that O — 0(1 — O) — 7(^1X2 +a;2a;3 -|-a;3a;i) 
and P = P [3(1 - O) - 7O], thus 
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PO - 30P 

: 7 [3PO(l - O) - 7PO2 - 3PO(l - O) 

+ 37P(a;ia:;2 + X2X3 + X3X1) 



-7^P 



O^ - 3(.xia;2 + X2X3 + X3X1) 



-^'^P Xl + a;| + ^3 - (a;iX2 + X2X3 + X3X1) 

(xi - X2)2 + (X2 - X3f- + (X3 - X2)^ 



< 0. 



So we need to find a function whose Lie derivative 



IS 7 



PO — 3OP . We notice that we can not in- 
tegrate it directly, however, we find out that the 
function 



<^ = 7 
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has the Lie derivative 

• PO - 30P 

(I> = 1- 
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< 0. 



(15) 



(16) 



As a result, we get a Lyapunov function. ^ = can 
happen at a point on the line xi = X2 = X3 or in 
the set (xi,X2,X3)|P = 0. 

B. Analysis on dynamics in the state space 

With the Lyapunov function constructed in the last 
subsection, we discuss dynamics in classified parameter 
space of the system here. 

1. 7 = 0: 

As ^ = only on the plane O = 1, all the tra- 
jectories will converge to this plane. Besides, on 
this plane, the value of the Lyapunov function 
(/) = 3 + InP will be a constant. Thus for each 
trajectory on the plane, the value of P will be a 



5 



constant. This means that the hmit set for any ini- 
tial point will be the intersection of the plane O = 1 
and the hyperboloid that P equals to a constant, 
which is a cycle on the plane. 

2. 7 < 0: 

As (j) = 7^ is non-positive now, the minimum 
value of (j) will not be zero if its initial value is 
not. Thus, in order to minimize the value of the 
Lyapunov function, all the trajectory will converge 
to one point on the line Xi = X2 = x^. So this case 
has a global stable equilibrium. 

3. 7 > 0: 

As = 7^ is non- negative now, the minimum 
value of (p will be zero. That is, all the trajectory 
will converge to the set {xi,X2,X3)\P = 0. In the 
neighborhood of P = 0, the terms of order X1X2, 
etc., in O asymptotically make a negligible contri- 
bution [3]. Thus O = 0{l-0) leads to O ^ 1 in 
the end. Finally, the limit set in this case is the set 
(a;i,a;2,a;3)|P = 0,O = 1. 

Three remarks are made here: 

• Our result is consistent with that in [3]. Further- 
more, wc give a full description on dynamics for the 
whole state space with the Lyapunov function. The 
energy landscape of the Lyapunov function on the 
plane xi + X2 + x^ (Fig. 2) gives a direct observa- 
tion on dynamics: (l)when 7 = 0, Fig. 2's case 

(1) shows the system has hamiltonian structure; 

(2) whcn 7 > 0, Fig. 2's case (2) shows the system 
has a global stable equilibrium; (3)when 7 < 0, 
Fig. 2's case (3) shows the limit set of the system 
is (a;i,a;2,a;3)|P = 0,0 = 1. 

• A similar analytic form of the Lyapunov function 
have been constructed when 7 > in [22, 23]. Com- 
pared with theirs, our construction is for the whole 
parameter space. Besides, we here provide a ex- 
plicit method to find this Lyapunov function. 

• In [24], C. W. Chi study the asymmetric May- 
Leonard system. Our construction method here 
may be able to be generalized to their system. 

V. CONCLUSION 

We have demonstrated that the Lyapunov function 
can be constructed in general 2-dimensional and two 
3-dimensional competitive Lotka-Volterra systems. For 
each example, we have shown dynamics in the whole state 
space with the Lyapunov function. The 2-dimensional 
case includes the bistable case and the model of May- 
Leonard has cycles as its limit set. Besides, in the 



Appendixes, wc have defined the generalized gradient 
system and discussed its coherence and generality with 
the classical gradient system. Furthermore, we notice 
that the construction method used in the model of May- 
Leonard may be able to be generalized to the asymmetric 
May-Leonard system in [24]. Thus the Lyapunov func- 
tion can be helpful to solve the limit cycle problems in 
general 3-dimensional case. 
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APPENDIXES 

In the Appendixes, wc first give the definition of the 
Lyapunov function and the construction framework in 
the Appendix I. Then we introduce the generalized gra- 
dient system in the Appendix II. Finally in the Appendix 
III, we give detailed calculation on all the dynamical 
parts in our framework of the third example. 

APPENDIX I. THE LYAPUNOV FUNCTION 

Definition 2. A smooth dynamical system is given by 

X = f (x) , (17) 

where x = (xi, • • • , a;„) with Xi, - ■ ■ , a;„ the n Cartesian 
coordinates of the state space, x = dx/dt and f : M" 
K". 

The conventional Lyapunov function for a given system 
(17) is defined as: 

Definition 3 (Conventional Lyapunov Function). [20] 

Let L : O ^ W be a function, where O is an open 
set in M" . L is a conventional Lyapunov function of the 
system (17) on O if 

• for a specified equilibrium x* in O, i(x*) = and 
L{x) > when x ^ x*; 

• i(x) = '^\^^0 for all X e O. 

La Salle has extended the conventional Lyapunov func- 
tion to include stable region by abandoning the positive 

definite requirement, but his generalization is too rough 
to lose stability information inside the stable region. 
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FIG. 2. The energy landscape of example (3) on the plane O = 1: (1) 7 = 0; (2) 7 = -0.1 < 0; (3) 7 = 0.1 > 0. 



Deflnition 4 (La Salle's Lyapunov Function). [25] 

Let L : O — !■ M be a function, where O is an open 
set in M". L is a La Salle's Lyapunov function of the 
system (17) on O if L{yi) = for all e O . 

Following the Lyapunov function used in [18, 19], here 
we give a more precise definition on the Lyapunov func- 
tion. 

Definition 5 (Lyapunov Function). Let (j> : M" i-> M &e a 
function, (p is a Lyapunov function of the system (17) 
if(f){x) = ^ |„ ^ /or all X e M" and (/)(x) only when 
X belongs to the union of the uj-limit sets UseRnCLi(s). 

In the following, we will briefly introduce our construc- 
tion framework. It is recently discovered during the study 
on the stability problem of a genetic switch [18, 26] and 
has been found wildly useful in biology [27]. The key 
result of the framework is a transformation from the n 
dimensional system (17) to the vector differential equa- 
tion (for simplicity, we only discuss the deterministic case 
in this paper, with noise strength being zero, general re- 
sults with randomness can be found in [18]): 



[5(x)+r(x)]i= -V0(x) 



(18) 



Here is a scalar function, the Lyapunov function. S 
is a semi-positive definite and symmetric matrix. T is a 
antisymmetric matrix. 

Symmetrically, if (5* + T) is nonsingular, the Eq. (18) 
can be rewritten as a reverse form 



X = - [D{^) + Q(x)] V0(x) , 



(19) 



where D is a semi-positive definite and symmetric matrix, 
Q is a antisymmetric matrix. 

From a physical point of view, S can be explained as 
a frictional force indicating dissipation of the potential 
energy, T as a Lorentz force and ^ as a potential of the 
system influenced by the two forces. Symbol D denotes 



the diffusion matrix indicating the random driving force, 
therefore for deterministic systems, D is free to choose. 
We make four remarks here: 



• As </> = i^Vc/) = -^^[5* -f- 
where r denotes transpose, 
be a Lyapunov function. 



T]± = -±^S± < 0, 
p in the Eq. (18) can 



• In this decomposition of the dynamical system, S 
can be considered as gradient part and T as rota- 
tional part. When S* = 0, it is a conserved sys- 
tem with first integral. If T is a scalar matrix at 
the same time, it is a Hamiltonian system where 
the trajectory would be a contour along the energy 
landscape of the Lyapunov function. When T ~ 0, 
it is a generalized gradient system defined in the 
next section. Thus both S and T can provide dy- 
namical information for a given system. 

• If the explicit form of the Lyapunov function is 
not obtained yet for a given system, such as the 
third example, we can solve the Eq. (19) by choos- 
ing a suitable form of D (or S) to make satisfy 
V X V(/) = and 4> < 0. Here Vx is a matrix gen- 
eralization of the "curl operator in 3 dimension" : 
(V X x)^^. = diXj — djXi. We have constructed the 
Lyapunov function by this method in the first two 
examples. 



• If one already has a Lyapunov function for a given 
system, we can obtain other dynamical parts [19]: 



f f 

V0 X f 

ff 



T = 



(20) 
(21) 



The corresponding explicit expression of the diffu- 
sion matrix D and the antisymmetric matrix Q can 
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be provided as well: 



As for the system (26) in Sec. Ill, the matrices are 



D = - 



f • f 



Q 



V0-f 

V</- X f 



/ + 



(V0 X f)' 



(Vet) ■ f ) {V(l> ■ Vet)) 



(22) 
(23) 



APPENDIX II. GENERALIZED GRADIENT 
SYSTEM 



P/xi 

S=\ a/x2 

P/X3 

xi/P 

D=\ x2/q 

X3//? 



C. Linear Cases 



T = 0, 



Q = 0. 



A. Definition 

Definition 6 (Generalized Gradient System). A gener- 
alized gradient system on is a dynamical system of 
the form 

X = -£>(x)V(/.(x) , (24) 

where ^ : M" i-^- R is a continuous differentiable scalar 
function and D{x) is a semi-positive definite and sym- 
metric matrix. 

By definition, when D is the product of a nonzero con- 
stant and the identity matrix, it degenerates to the clas- 
sical gradient system in [20]. The potential gradient of 
the system (24) is anisotropic, different from that of the 
gradient system. Such anisotropic system has been ob- 
served in real systems, like the Fourier's equation in [28] . 
T = is equivalent to Q = as 

r(x) = 1 [(D(x) + Q(x))-i - (I)(x) - Q(x))-i] . 



B. Matrices S, T, D and Q of the first two examples 

For the given system (2) in Sec. II, it is not a gradient 

system as the curl of the vector field V x x = axi — f3x2 ^ 
0. But we can calculate the matrices using the results 
obtained in Sec. II. A: 

D being semi-positive definite and symmetric and T = 
indicate that the system (2) is a generalized gradient sys- 
tem with zero rotational part, and trajectory will not con- 
tour along the energy landscape of the Lyapunov func- 
tion. Besides, S is singular only on the coordinate axis in 
this system. It means the dissipation is infinite and thus 
the trajectory will stay on the axis once reaching it and 
approach the equilibrium E^q = (&i,0) or Eq^ = (0,62)- 



A linear autonomous dynamical system is given by the 
following ordinary differential equations: 

X = Fx , (25) 

where x = (xi, • • • , .x„) with xi, ■ ■ ■ , a;„ the n Cartesian 
coordinates of the state space, x = dx/dt and F is a con- 
stant matrix. To ensure the independence of all the state 
variables, we require the determinant of the F matrix to 
be finite: det{F) ^ 0. 

To illustrate the coherence and generality of the gen- 
eralized gradient system in the linear cases, we first 
mention that a linear system (25) is a gradient system 
X = — if and only if its F matrix is symmetric: 

• A gradient system x = — V<ji has d(f)/dxi = 
~T,"^-^F^jXj, then V x Vcj) = leads to the F ma- 
trix being symmetric: 

• If a linear system (25) has a symmetric F matrix, 
then by setting d(j)/dxi = —'EJ^iFijXj, the solution 
of ^ exists, and we can rewrite (25) as x = — V^. 

But a linear system (25) can be a generalized gradient 
system when the F matrix is asymmetric. Such systems 
have nonzero curl, that is V x x 7^ 0. We give an example 
of 2 dimensional linear generalized gradient system in the 
following. 

Example 4. This example is given by [20]: 

We set a Lyapunov function to be ^ = — X\X2 as its 
Lie derivative 

(j) = -3a;i - {xi - 2x2f < . 
Then the system (26) can be rewritten as 




Therefore, the original system (26) is a generalized gra- 
dient system by definition, but not a gradient system as 
its F matrix is asymmetric. 
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III. MATRICES S AND T FOR THE MODEL OF 
MAY-LEONARD 



We notice that 1 — q; = — (1 — /3) = in the case 
of 7 = 0. Thus we have 



In this section, we calculate S and T by the Eq. (20) 
for the model of May-Leonard system. In Sec. V, we 
have obtained the Lyapunov function: 

1. When 7 = 0: 

(f> = 30-lnP 
= 3{xi +X2+ xs) - ln(a;ia;2a;3). 



2. When 7 7^ 0: 



(f> = l 



= 7 



03 



X1X2X3 



{Xi + X2+ X3) 



3 ■ 



Thus we do calculation separately for this two cases in 
the following. 

1. When 7 = 0: 
As 



3x1 



3x2 — 1 

dx2 X2 

d<j> _ 3x3-1 

9X3 



(28) 



f f 

3[1 - {Xi+X2+Xs)]^ 



X -I I I Xo 



-I . (29) 

-r ~<~ -^3 

Notice on the plane Xi + X2 + X3 = 1, S' is zero 
matrix, thus on the plane the system is conserved 
when 7 = 0. 

As for T: 

V(t>xf 



f f 

'3xi-l, 



3x, —1 



3x3 



(30) 



Since T is antisymmetric, we just need to calculate 
the elements T12, T13 and T23 of the above 3x3 
matrix. As we have proved all the trajectory will 
converge to the plane X1+X2+X3 = 1, we calculate 
the elements on the plane below. We first calculate 
T12. 



nn 3a;i - 1 . 

T12 = — X2 

X\ 

2>x\ — 1 



3a;2 - 1 . 

x\ 

X2 



(1 - PX\ — X2- aX3)X2 
Xi 

3X2 - 1 



X2 



-(1 - Xi- O.X2 - PX3)XI . 



(31) 



Tl2 = 



H [ix2 - Xi) + {X2 - X3)]{X3 - X2) 

X2 



— [{xi - X2) + (xi -a;3)](a;3 - xi) 

X\ 



We calculate T13 and T23 similarly. 



3a;i - 1 . 3a;3 - 1 . 
Ti3 = 0:3 Xx 



(32) 



3a;i - 1 



X3 



(1 - aX\ - 13X2 - X3)X3 
Xl 

3X3-1, 



X3 



-{1 — Xl — ax2 — Px3)xi 



l3-a 



Xl I 



— [{Xl - X2) + {xi - X3)]{xi - X2) 
Xl 



H [{^3 - a;i) + (a;3 - 0:2)] (0:2 - 0:3) 

X3 



^ 3X2-1. 3X3-1. 

T23 = 2:3 X2 



(33) 



X2 
3X2 - 1 



X3 



(1 - aXi - (3X2 - X3)X3 
3X3 - 1 , 



X3 



-(1 - f3xi - X2- aa;3)a;2 



X2, 



— [{X2 - Xl) + {X2 - X3)]{xi - X2) 
X2 



H [{^3 - X-i) + {X3 - X2)]{X3 - Xl) 

X3 



(34) 



Thus we calculate out each element of matrix T on 
the plane X1+X2 + X3 = 1. 

2. When 7^0: 
As 

= l-§i[x2X3{xi +X2 + X3) - 3.xi.T2a;3] 

§^=l-§s[xiX3{xi+X2+X3) -3X1X2X3] , (35) 

Ml = l-as[xiX2{xi +X2+ X3) - 3a;ia;2a;3] 



5 = -^7 



f f 



20* 



{Xi - X2f + (X2 - X3f + {X3 - X2f 



±1+ X2 



^2 



(36) 



Since on the plane xi + X2 + X3 = 1, S 
is not zero matrix except on the limit set 
{xi,x2,x3)\P = 0, = 1. Thus on the plane the 
system is dissipative when 7 0. 
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As for T: 



T 



V0 X f 
f f 



7cF 



3x3 



(37) 



Again, as all the trajectory will converge to the 
plane O = xi + X2 + xy, — 1, we just calculate the 
elements of above 3x3 matrix on the plane: 



1 3x^ . 1 3x^ 



(38) 



Here we use Tij to denote the matrix elements so 
that they can be distinguished with the matrix ele- 
ments Tij in the case of 7 = 0. Since l — a= 
and 1 — /3 = — in the case of 7 7^ 0, we get 

712, and T23 with similar calculation: 



X2 ' 



[{X2 - Xi) + {X2 - X3)][{1 - a)x2 + (1 - I3)X3] 



X2 
Xi 



[{xi - X2) + {xi - X3)][{1 ~ a)x3 + (1 - /3)xi] 

(39) 



We calculate T13 and T23 similarly. 



X3 ' 



[{X3 - Xi) + {X3 - X2)][il - a)x2 + (1 - I3)X3] 



X3 

Xi 



[{xi - X2) + (xi - X3)][{1 - a)xi + (1 - I3)x2] 

(40) 



X2 
X3 ' 



[{x3 - xi) + {x3 - X2WI - a)x3 + (1 - /3)xi] 



X3 

X2 ' 



[{x2 - xi) + {x2 - X3)][{1 - a)xi + (1 - I3)X2] 

(41) 



Thus we calculate out each element of matrix T on 
the plane xi + 2:2 + 2:3 = 1. 
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